Boundary conditions for quasiclassical equations in the theory of superconductivity 
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In this paper we derive effective boundary conditions connecting the quasiclassical Green's function 
through tunnel barriers in superconducting - normal hybrid (S-N or S-S') structures in the dirty limit. 
Our work extends previous treatments confined to the small transparency limit. This is achieved by 
an expansion in the small parameter r _1 = T/2(l — T) where T is the transparency of the barrier. 
We calculate the next term in the r _1 expansion for both the normal and the superconducting case. 
In both cases this involves the solution of an integral equation, which we obtain numerically. While 
in the normal case our treatment only leads to a quantitative change in the barrier resistance Rb, 
in the superconductor case, qualitative different boundary conditions are derived. To illustrate the 
physical consequences of the modified boundary conditions, we calculate the Josephson current and 
show that the next term in the r _1 expansion gives rise to the second harmonic in the current-phase 
relation. 
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I. INTRODUCTION 

In recent years there has been a great deal of work, 
both experimentally and theoretically, on the transport 
properties of mesoscopic superconducting - normal (N- 
S) hybrid structures. Various new effects have been ob- 
served including zero-bias anomalies in the differential 
conductance Q], peculiar dependence on the magnetic 
field, re-entrant temperature behaviour All these 
phenomena are caused by the interplay of Andreev scat- 
tering at N-S interfaces and phase coherence in the meso- 
scopic region and are manifestations of the proximity ef- 
fect, whereby superconducting correlations are induced 
in the normal region. 

Among the various theoretical techniques used to deal 
with the above effects, the quasiclassical method has re- 
vealed itself as one of the most powerful approaches (see, 
for instance, [[§ - @)- This technique enables the study of 
thermodynamical and kinetic properties of superconduc- 
tors (see H and references therein), whose dimensions 
significantly exceed the Fermi wave length \p = 2n/kp 
and where quantum size effects can be neglected. In the 
case of hybrid structures the quasiclassical equations of 
motion must be supplemented by appropriate boundary 
conditions in order to match the quasiclassical Green's 
function g at N-S interfaces. These have been derived by 
Zaitsev Q in a general form which is valid in both the 
clean and dirty limits and take the form |10[] 
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Here a and s are the antisymmetric and symmetric 
parts of the supermatrix Green's function g, i.e. 



g(±n) = s±a, 



(2) 



where /i = cos(8) = p z /pf and 9 is the angle between 
the velocity p/m of an incident electron and the vector 



normal to the interface. The reflection (R) and transmis- 
sion (T) coefficients depend on \i and are connected via 
the unitarity relation 



R((i)+T(n) = l. 



(3) 



As shown in Q , the antisymmetric part a is continuous 
at the interface, while the symmetric part s experiences 
a jump determined by the commutator on the right hand 
side of eq.(Q) 



[S2,Sl] = S 2 Si - SiS 2 - 



(4) 



The matrix elements of g are the retarded (advanced) 
Green's functions g R ( A ) and the Keldysh Green's func- 
tion g 
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Following the usual convention, we denote two-by-two 
matrices in the Nambu space by a "hat" (<?) symbol and 
the four- by- four supermatrices by a "check" (g) symbol. 
The Keldysh Green's function g describes the kinetic ef- 
fects and, by exploiting the normalization property, 



99 = I 

is related to the matrix distribution function / via 

9 = 9 R f-fg A - (6) 

In what follows we will adopt the convention intro- 
duced by Larkin and Ovchinnikov, according to which 
the matrix / can be chosen to be diagonal, with / = 

/ol + fz&z- 

The boundary conditions (Q) are valid in both the clean 
and dirty limit. We recall that in the dirty limit (I <C 
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£n,Sj I is the mean free path, £jv,,s are the coherence 
lengths in the N and S regions, respectively) the angular 
dependence of the matrix g can be taken into account 
by keeping only the first two terms in an expansion in 
Legendre polynomials P n (fi). The first two terms in the 
expansion of g can be then related to each other. Far 
away from the interface this relation can be written down 
as | 



The layout of our paper is as follows. As a prelimi- 
nary step, in the next section, we consider the case of 
two normal regions separated by a barrier. In the follow- 
ing section we turn our attention to the superconducting 
case. In section IV we compute the Josephson current in 
the presence of the modified boundary conditions. Our 
conclusions are finally stated in section V. To simplify 
the notation in the following we will drop the explicit fi 
dependence of r, and only reinstate it whenever neces- 
sary. 
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1. Schematic picture of the structure studied in the 



The subscript oo means that \z\ 3> I (the interface is 
located at z = 0, see Fig.l). However, close to the inter- 
face, one should keep all the terms in the expansion of g 
because the coefficients R and T entering the boundary 
condition (Q) depend on /i. In the dirty limit, all higher 
terms ( n > 2) in the expansion of g decay exponentially 
with z/l, and it would be desirable to obtain a matching 
condition at the interface which involves the asymptotic 
functions cioo and only. Such a problem was con- 
sidered by Kupriyanov and Lukichev [|l0|. They showed 
that the boundary condition (pi) reduces then to 



«2e 



(8) 



II. NORMAL CONDUCTOR WITH A BARRIER 

The case of two normal conductors (N-I-N) separated 
by a barrier was analyzed by Laikhtman and Luryi [ pd| , 
p"2| . We present here a derivation of the effective bound- 
ary conditions which differs from that presented in pd| ] 
and is applicable also to the more general case when one 
or both conductors are in the superconducting state (S- 
I-N and S-I-S). In the case of normal conductors, eq.(Q) 
implies the following boundary condition 



/ 2 (/i)=T/ 1 ( M )+ J R/ 2 (- / x) 



(11) 



where fi.2 = (fa + /z)i,2 is the usual distribution func- 
tion, as it is clear by recalling the definitions of the func- 
tions /o and f z in eq.(0) (see also 0). Note also that, in 
the normal case, eq.(O) can be easily derived by means 
of simple counting arguments. By considering the sym- 
metric (s) and antisymmetric (a) parts of /, one rewrites 
the boundary condition (pT[) in the form 



ro(0) = [s] = « 2 (0) - si(0), 



(12) 



where we have made use of the continuity of the anti- 
symmetric part a at the boundary. 

In order to find a spatial dependence of the distribution 
function, one has to solve the Boltzmann kinetic equa- 
tion, in the relaxation time approximation in terms of s 
and a reads 



-a/n — —b 



(13) 



< fj,/r(fx) > [s 2oo ,si 



(9) 



Here r(/i) = 2R(^)/T(fi) and the angular brackets 
mean the angle averaging 



< p/r(ji) >- 



d\i n/r{p). 



(10) 



The main aim of the present paper is to show that 
eq.(|^) is valid only to lowest order in an expansion in 
the small parameter r~F i.e. in the limit of a low trans- 
parency barrier, and consequently in the case of an arbi- 
trary barrier transmission, eq.(|^) can be used only for a 
qualitative description. 



H 2 b' =< s > -s 



(14) 



where to make our notation more compact s' = ld z s 
and, we assume for simplicity that the mean free path, 
I, is the same on both sides of the barrier. Eliminating s 
from eqs.( |l3||L4] ), we can write the equation for b as 

(15) 



H 2 b" -b=- <b> +B S{z/l) 



which is valid for all z. The function Bo(fi) = Bq 
determines the jump in the derivative b' at z = 

B Q = n 2 [b'\ = [< s > -s] = rM&(0,/z)- < rfib(0,fi) > 

(16) 
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where the last equality follows from eq. ( |l2| ) . We single 
out the asymptotic part of b(z) and s{z) 



b{z) — boo + 5b(z), s(z) = s ao + 5s(z) 



(17) 



in such a way that the functions boo and Soo do not 
depend on \i and are connected to each other via the 
relation (see eq.(|l3|)) 



-s'(±oo). 



(18) 



It follows from eq.(|l4]) that the average of < /i 2 b > 
does not depend on the coordinate (this amounts to the 
conservation of the current) and therefore can be evalu- 
ated from the asymptotic part 



< [i 2 b >= ^ =< M 2 &(0, M) > 



(19) 



As a result eq. ([L5|) can now be written for the function 
5b decaying at \z\ — > oo 



H 2 Sb" -5b = - <5b> +B 5(z/l). 



(20) 



By performing a Fourier transform, we obtain for the 
Fourier component 6b q 



Sb q = m q (< Sb q > -B Q ) 
and hence for the average 



(21) 



By a little manipulation of the boundary condition (|1 
we get 



rnb(0,n) = [soo] + [5s] 



(25) 



where [soo] = s(oo) — s(— oo). 

Here the jump [5s] = <Js2(0) — <5si(0) and taking into 
account the continuity of 5b at z = 0, we obtain from 



[5s] = f 

J — ; 



dzSb(z) = 5b, 



(26) 



where qo = 0. By using eq .(pi| ) to express 5b qo and 
taking into account both eqs.(|25|) and (|26]), after angle 
averaging, we arrive at the desired effective boundary 
condition 



3 < r(^i)^ 3 b(fi) > boo = [s c 



(27) 



Once the solution of the integral equation b(fi) is 
known, eq.(|27]) provides a relation between and the 
jump of the symmetric part [soo]- Such a relation is use- 
ful in evaluating the current density j, which can be ex- 
pressed in terms of b^ as 
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a 

2d 



de b r 



(28) 



< 5b q >- 



1— < m q > 



< m q ii B Q >, 



(22) 



where m„ 



1/(1 + pi 2 q 2 ), < m q >= arctan(q)/q. 
By performing the inverse Fourier transform, we find the 
magnitude of b{z, ji) at z = 



where a = 2e 2 N vpl/3 is the conductivity, No being 
the single particle density of states per spin. If the sym- 
metric part of the distribution function is in equilibrium 
on both sides of the barrier, then 



[s c 



th[(e + eV)0\-th[(e-eV)0] 



(29) 



1— < m q > 



b(0, fi) = boo + I TT-mq 



< m q r(fi)fi 3 b(0, n) > -r((i)fib(0,(ji) 



(23) 

Eq.(p3|) is an integral equation for b(0,fi). It can be 
rewritten in the form 



b(/x)(l +r(ju)/2) = 1 + / diMir(nx)niK,(n,ni)b(ni) 
Jo 

(24) 



where 2V is the voltage drop across the barrier, (3 — 
1/2T. As a result, the barrier resistance per unit area 
becomes 



3/ 

R b = (2V/j) = -< r(n)n a b(n) > 
a 



(30) 



This formula gives a relationship between the barrier 
resistance Rb and the reflection and transmission coeffi- 
cients. The function 6(/x) must be found from eq. (|24]) . In 
the 3-dimensional case eq.(p0|) reads 



R b = R bQ < r{n)n z b{ii) > 



(31) 



where 



where 



6(a0=6(0,m)/& o 



dq 



2 2 

1 Mi 



RbO = 



9ir 2 h 



The integral equation ( p4[ ) must be solved numerically, 
but before describing the numerical solution, it instruc- 
tive to consider few limiting cases, where an analytical 
approach is possible. 
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A. Weak barrier 



Mm) = x(m)/(?-aO- 



(39) 



This means that r€l (this condition should be ful- 
filled for angles not too close to 7r/2). Then, it is seen 
from eq.(E4) that 



and 



Hp) « 1 



Rbw = RbO < r/i > 



(32) 



(33) 



For example, in the case of a thin barrier, modelled by 
a delta-like potential at z = 0, we have that 



and we obtain 



r(/0 = S 2 /M 2 



(34) 



(35) 



where s = \/2UbW / '(v fU)\ U b , w, and «f being the 
barrier height, barrier width, and Fermi velocity, respec- 
tively. 



B. Strong barrier 



This means that r«l. To lowest order the solution 



to eq.(24) is obtained as 



(36) 



where C is a /i-independent constant which can be 
determined by exploiting the fact that the average < 
/i 2 6(0,/i) > does not depend on z. In fact, as follows 
from eq. (|l9|) 



Substituting (|39| ) into eq.(g4|), we obtain an integral 
equation for x(p) 

3^r < M /r> ~ 1 = "^ + i ^i/C(M,Mi)x(Mi). 

(40) 



It follows from eq.(40) that the function x is of the or- 
der 1, i.e. it is small compared to C in eq. (|36|) . Therefore, 
in order to obtain the effective boundary condition con- 
necting boa and [soo] m the general case, we must solve 
the integral equation (p4]). 



C. Numerical results 

In the appendix we discuss how the integral equation 
is solved and here we merely illustrate the numerical re- 
sults. Figure 2 shows the behaviour of the function b(fi) 
for different values of the parameter s. We see that at 
small and large s the numerics yield the expected be- 
haviour discussed above. It is perhaps worth noting that 
for a planar barrier the Landauer formula yields 



Rbarrier — ^ n 



h ( 2?r 



2e 2 \k F 2 



< 



>" 



(41) 



which agrees with eq.(p7|) and is smaller than eq.(|33|) 
by a factor of 8/9. Figure 3 shows the behaviour of the 
barrier resistance as a function of s normalized to the 
s = resistance. 

In the next section, we now turn to the analysis of the 
more general case of a S-I-N or S-I-S interface. 



C= {3<Li/r>y 1 
which lead to a barrier resistance equal to 

R bs = R b0 {9< fJL/r>y 1 . (37) 



The above results coincides with that obtained by Ref. 
Of, as can be noted from eq.(0). For a thick barrier, we 



get 



Rbs = RbO^S 



(38) 



We note that if we were to use, wrongly, eq. ( |37j ) in the 
case of a weak barrier, we would obtain the expression 
( |38| ) instead of (|35|). The ratio of these two results, 9/8, 
is close to 1. However, this ratio does not contain any 
small parameter. This conclusion is in agreement with 
that of Ref. [O. 



One can obtain a correction 6i(^t) to the solution (36) 
To see this, we seek a solution in the form 




FIG. 2. The solution of the integral equation, b(fi) is 
plotted as function of [t, for different values of s: s — solid 
line, s — 1 dashed line, s = 10 long-dashed line. The mesh 
size in space is 1024 points. 
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FIG. 3. Barrier resistance Rk versus s. 



III. GENERAL CASE: S-I-N AND S-I-S 

In this section we analyze boundary conditions for the 
quasiclassical matrix Green's function g in the dirty limit. 
We will start with the boundary conditions for g obtained 
by Zaitsev in the general case Jj| and, as in the preceding 
section, find a relation connecting the asymptotic values 
of g at \z\ ^> I on the two sides of the barrier. In contrast 
to the previous case, the presence of anomalous compo- 
nents in the matrix Green's function leads to a nonlinear 
equation governing g. In order to simplify the problem, 
we note that in the case of a weak barrier (i.e. R — > 0, 
T — > 1) the boundary condition for the symmetric parts 
sj.,2 reduces to a continuity condition at the interface. It 
is then physically relevant to confine our analysis to the 
most interesting case of a strong barrier and obtain a re- 
lation between a and s using the expansion in the small 
parameter r _1 . 

By assuming that r ^> 1, eq.(||) can be cast in the form 



1 



& (0>M) ~ o — [«2,Sl] 
zr/i 



1 , 

1 - 7T 
2r 



(42) 



where as in the previous section we introduced the 
function b = a/fx. Eq.([42"|) is valid up to terms of the 
order r~ 2 . From the equations for s and b (see below), 
one can obtain the relation corresponding to eq. ( fL9| ) , with 
b replaced by the supermatrix b. We then proceed as be- 
fore by writing s and b as the sum of a fast decaying and 
an asymptotic part 



b = b n 



5b 



(43) 



where we assume 5s <C Soo- By multiplying eq.([42|) 
by /i 2 and using the representation (J43|), we perform the 
angle average of eq. pa) to obtain 



g&oo = [s 2(x5) S loo ] ^< — > - < ^ > (S 2c 



2r 



Sic 



< ir ([s 2 oo,^Si] + [<5s 2 ,Sloo]) > 

2r 



(44) 



Here 5§2 i = 5s2,i(Q ± )- The problem is then reduced 
to the calculation of the functions 5§2,i- We start by 
writing the equation for g in the space interval < \z\ <€. 
£,n,s- It is then sufficient to retain only the gradient 
term and the collision integral in the self-consistent Born 
approximation for the impurity scattering. As a result 
the equation for g reads (|] 



2 W?' = 9 <<?> - <3> 9- 



(45) 



By rewriting eq. ( (45| ) in terms of b and s we get 

2fi 2 b' = s<s>~<s>s (46) 



2s' = b<§> — <3>b 



(47) 



together with the conditions deriving from the normal- 
ization condition gg = 1, 



1, sb + bs = 0. 



(48) 



Using the expansion (43) we obtain the equations for 
the deviations 5b and 5s in the form 



(49) 



(50) 



(51) 



(52) 



which follow from (|4§|). From eqs. ( |49| - |50| ) we finally get 
the equation for 5b 



li 2 5b' — —Soo (5s— < 5s >) 

5s' = Soodb s'ac = -Soo&oo. 

where we have used the relations 

Sbsoo + SooSb = 0, Soo5s + Sssoo = 

Soo Soo — 1 



H z 6b" -5b= - <5b> +B 5{z/l). 



(53) 



One can easily check that the matrix Bq is connected 
to the Fourier component 5b qo (where qo = 0) by 



B =< 5b q , > -5b„ 



From eq.(53) we find the Fourier components 

„2 



5b q = Tflq 



1— < m q > 



< m q Li 2 B Q > -B Q 



(54) 



(55) 



and the value of 5b(0, fi) at z = reads 
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Sb(p,n) = h(p,/i)-h 



dq 



2n V 1- < m q > 



< m q fi 2 Sb qo > —5b, 



(56) 



To close the above equation, we need to connect &(0, pi) 
and 600. To lowest order in r _1 , the boundary condition 
© yields 



6(0, /i) » (2r/i)- 1 [s 2oo , 



•Sic 



&oo ~ 3 < /i/2r > [s 2oo , sioo] ■ (57) 
which when substituted into (p6[) yields the equation 



current. To this end we rewrite eq.(|6l]) in the following 
way 



1 



b = A [92,9i) + B[g 2 ,gi]{g2,gi} 



(62) 



where we have identified the symmetric part s with 
the Green's function g and dropped the oo suffix and the 
curly brackets indicate the anticommutator 



{fa,, 9i } = 929i + fafa- 
The constants A and B can be read off from eq.(|6l|) 



A =< fx/r > -2 < fi/4r 2 > 



and 



1 



3/ir < fi/r > 



-1 = 



2n \ 1— < m q > 



If we seek the solution in the form 



Sb a 



(58) 



(59) 



then the function \ satisfies the integral equation fiO) 
and the Fourier component Sb qo is related to <5si,2 by 
integrating eq.(p0|) from — oo to and from to oo, to 
yield 



Sh = s 2o o8b qo /2, Ssi = -siooSb qo /2. 



(60) 



Substituting eqs.(]60j), (|9|), and (|5j) into eq.(|4j) we 
finally obtain 



r- - i ( M A* /- 



o MX M ^ r- - - - 1 

O < — >< — > Sloo, S2ooSlc5oS2oo 

2r 2r 



(61) 



The above equation is the effective boundary condition 
for the matrix g in the dirty limit. The first term in ( |6l| ) 
coincides with the boundary condition obtained in Rcf. 



B =< fi/Ar 2 > -3 < /i%/2r >< /j,/2r > . 

In deriving eq.(|62]) we have made use of the normaliza- 
tion condition gg = 1. The current through the junction 
is determined by the formula 



I = - 



I6el 



deTr(a z b) 



(63) 



where b is the appropriate Keldysh component. In the 
absence of a voltage across the junction, the Keldysh 
component of the supermatrix g reduces to 



(64) 



where /o = 2tanh(e/2T) is the equilibrium distribution 
function. The Keldysh component of the product of the 
commutator and the anticommutator reads 

{[fa,9i]{fa,9i})k = [g2>9i]{fa>gi}k + [fa,gi}k{g£,9f} 

(65) 

with the Keldysh component of the commutator 

[fa,gi]k = fo([g 2 H ,g^}-[g 2 A ,gt}) (66) 

and of the anticommutator 

{fa,9i} k = hiig^gf} - {gtg?})- (67) 

To calculate the current we represent the matrices in 
the Nambu space as 



IV. THE JOSEPHSON EFFECT 

As a simple application of the above boundary con- 
dition, we now derive an expression for the Josephson 



§ R(A) = G R(A) &z + lF R(A) {cos{(j)i2)ay + sm{(j)12)dx) 

where 4>i^ are the phases of the superconducting order 
parameter on the two sides of the junction. The current 
in eq.(|63|) can be written then as the sum of three terms 
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I A = -Uo^sinifa - 4> 2 ) / def Q {F*F* - F^), 

J — OO 

(68) 



= -Uo,B2sin((f>i - 4> 2 ) 

/OO 
def (F*F*G*G* - F A F 2 A G?G£), (69) 
-OO 

and 

(2) 

I B = -Ho^Bsin(2(4> 1 - <j) 2 )) 

/OO 
A/o((F«F«) 2 - (ifi^) 2 )- (70) 
-OO 

In the above formulae Iq^a = (sNqVf /2)A, Jq,b = 
(eNoVF/2)B, e is the electron charge, -ZVg the single par- 
ticle density of states per spin, and vf the Fermi velocity. 
By using the expression for G R ^ and F R ^ at equilib- 
rium 

G R(A) = e F R(A) = A 

v/(^±) 2 -A2' V(e±) 2 -A 2 ' 

where e± = e ± i0 + , and assuming that the gap A is 
equal on both sides of the junction, we obtain, at T = 0, 
the following result for the current 

I = eNoVpAn [(2A + B)sin(<j)) - Bsin(2<j))} (71) 

where <fi — <\> 2 — 4>\ . Note that by confining ourselves to 
the lowest order in r _1 , we would obtain for the Joseph- 
son current the standard result of tunneling theory 

/ = (e 2 N v F < fi/r >)^sin{0) = R^^aintf) (72) 

with Rb s given by eq.(|38"|). Allowance for higher order 
terms in the barrier transparency leads then to higher 
harmonics in the current phase relation. This result has 
a simple physical interpretation. We know that in the 
case of a superconductor - normal metal - superconduc- 
tor structure, the Josephson effect manifests itself with 
a triangular shape of the current - phase relation. For 
this case the Fourier decomposition has an infinite num- 
ber of harmonics. Hence it is clear that higher order 
terms in the r _1 expansion must possess harmonics of 
higher order. It may be useful at this point to notice 
that very recently Josephson current measurements have 
been carried out in structures in which it is possible to 
control experimentally the barrier transparency. The ex- 
perimental results show indeed that by tuning the barrier 
strength it is possible to observe higher harmonics |l3) . 



V. CONCLUSIONS 

We have analyzed the boundary condition for the qua- 
siclassical matrix Green's function g at the S-N or S-S' 
boundaries in the dirty limit. Effective boundary con- 
ditions have been derived with the aid of an expan- 
sion in the barrier transmittance. The first term of 
the expansion reproduces the results previously derived 
by Kupriyanov and Lukichev pJ. In the normal case, 
the boundary conditions for a barrier of arbitrary trans- 
parency may be obtained from the solution of an integral 
equation. In the superconducting case, due to the non- 
linear nature of the equations, we have been able to com- 
pute the next-to-leading term in the r _1 expansion. In 
this case the evaluation of this term also entails solving an 
integral equation. To illustrate the physical consequences 
of the modified boundary condition, we have calculated 
the Josephson current for a tunnel junction between two 
superconductors at equilibrium and shown that higher 
order harmonics arise. Given the relevance that the low- 
est order boundary expression (cf. eq.(^)) has played in 
the study of mesoscopic normal-superconducting hybrid 
structures, our result calls for a reexamination of the var- 
ious effects occurring in mesoscopic hybrid structures. In 
particular, one may envisage a straightforward general- 
ization of the circuit theory of Nazarov (m) to allow for 
the new boundary condition. This will be the subject of 
a future investigation. 
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APPENDIX A: SOLVING THE INTEGRAL 
EQUATION: TECHNICAL DETAILS 

The variable fi is discretized \n — iA^ with i = 1, ...,N 
and A p = 1/N. By introducing a vector bi = b(fii), the 
integral equation ( p4| ) acquires the matrix form 

JV 

6,(1+^/2) = 1 + Y t K ij r j n j b j (Al) 
i=i 

where = r(/ii) and 

/oo f^ 2 ^/ 2 
— m q (ij,i)m q (fj,j)- 3 ■ (A2) 
- oo 27T I- <m q > 

The solution bi can then be obtained as 

N 

bi = J2( A ~'h (A3) 
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where the matrix Ay is given by 

Ay = (1 + n/2)S lj - A^KijrjfXj. 



(A4) 
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